Setup and Visualize
Load the rwrfhydro package.
library("rwrfhydro")
## To check rwrfhydro updates run: CheckForUpdates()
# set your own path
rlFile <- '~/WRF_Hydro/DOMAIN_library/Boulder_Creek_100m_1km_2sqkm_full_2015_09_03//Route_Link.nc'
# currently there is no "gages" variable
ncdump(rlFile)
## File: ~/WRF_Hydro/DOMAIN_library/Boulder_Creek_100m_1km_2sqkm_full_2015_09_03//Route_Link.nc
## ( NC_FORMAT_NETCDF4 ):
## dimensions (1):
## linkDim = 403 ;
## variables (24):
## int link(linkDim) ;
## link:long_name = "Link ID" ;
## int from(linkDim) ;
## from:long_name = "From Link ID" ;
## from:coordinates = "time lat lon" ;
## int to(linkDim) ;
## to:long_name = "To Link ID" ;
## to:coordinates = "time lat lon" ;
## float lon(linkDim) ;
## lon:long_name = "longitude of the start node" ;
## lon:units = "degrees_east" ;
## lon:standard_name = "longitude" ;
## float lat(linkDim) ;
## lat:long_name = "latitude of the start node" ;
## lat:units = "degrees_north" ;
## lat:standard_name = "latitude" ;
## float alt(linkDim) ;
## alt:long_name = "Elevation in meters at start node" ;
## alt:standard_name = "height" ;
## alt:units = "m" ;
## alt:positive = "up" ;
## alt:axis = "Z" ;
## int type(linkDim) ;
## type:long_name = "Link type" ;
## type:coordinates = "time lat lon" ;
## int order(linkDim) ;
## order:long_name = "Stream order (Strahler)" ;
## order:coordinates = "time lat lon" ;
## float Qi(linkDim) ;
## Qi:long_name = "Initial flow in link (CMS)" ;
## Qi:coordinates = "time lat lon" ;
## float MusK(linkDim) ;
## MusK:long_name = "Muskingum routing time (s)" ;
## MusK:coordinates = "time lat lon" ;
## float MusX(linkDim) ;
## MusX:long_name = "Muskingum weighting coefficient" ;
## MusX:coordinates = "time lat lon" ;
## float Length(linkDim) ;
## Length:long_name = "Stream length (m)" ;
## Length:coordinates = "time lat lon" ;
## float n(linkDim) ;
## n:long_name = "Manning's roughness" ;
## n:coordinates = "time lat lon" ;
## float So(linkDim) ;
## So:long_name = "Slope (%; drop/length)" ;
## So:coordinates = "time lat lon" ;
## float ChSlp(linkDim) ;
## ChSlp:long_name = "Channel side slope (%; drop/length)" ;
## ChSlp:coordinates = "time lat lon" ;
## float BtmWdth(linkDim) ;
## BtmWdth:long_name = "Bottom width of channel" ;
## BtmWdth:coordinates = "time lat lon" ;
## float LkHZArea(linkDim) ;
## LkHZArea:long_name = "Maximum lake area (sq. m)" ;
## LkHZArea:coordinates = "time lat lon" ;
## float LkMxH(linkDim) ;
## LkMxH:long_name = "Maximum lake elevation (m)" ;
## LkMxH:coordinates = "time lat lon" ;
## float WeirC(linkDim) ;
## WeirC:long_name = "Weir coefficient" ;
## WeirC:coordinates = "time lat lon" ;
## float WeirL(linkDim) ;
## WeirL:long_name = "Weir length (m)" ;
## WeirL:coordinates = "time lat lon" ;
## float OrificeC(linkDim) ;
## OrificeC:long_name = "Orifice coefficient" ;
## OrificeC:coordinates = "time lat lon" ;
## float OrificeA(linkDim) ;
## OrificeA:long_name = "Orifice croess-sectional area (sq. m)" ;
## OrificeA:coordinates = "time lat lon" ;
## float OrificeE(linkDim) ;
## OrificeE:long_name = "Orifice elevation (m)" ;
## OrificeE:coordinates = "time lat lon" ;
## float time() ;
## time:standard_name = "time" ;
## time:long_name = "time of measurement" ;
## time:units = "days since 2000-01-01 00:00:00" ;
##
## // global attributes (2):
## :featureType = "point"
## :history = "Created Thu Sep 03 09:30:05 2015"
# This function returns a function which plots the data
PlotRouteLink <- VisualizeRouteLink(rlFile)
PlotRouteLink() # basic plot
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=39.940325,-105.430153&zoom=8&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
args(PlotRouteLink) # what are the options? unfortunately I havent really documented these.
## function (location = c(lon = mean(range(rl$lon)), lat = mean(range(rl$lat))),
## indices = FALSE, comIds = FALSE, zoom = if (length(gageZoom)) 13 else 8,
## source = "google", maptype = "terrain", padPlot = 0.01, linkColor = "blue",
## textColor = "darkred", gageColor = "cyan", gageShape = 9,
## gageZoom = NULL, plotPath = NULL, plotType = "pdf", doPlot = TRUE)
## NULL
PlotRouteLink(zoom=10, comIds=TRUE) # links on the edge of the domain arent fully plotted... because there's no "to" lat/lon.
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=39.940325,-105.430153&zoom=10&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
Create the gages variable
Define the synthetic gages to be used for the domain. In this example put gages on ALL the links.
rlLink <- ncdump(rlFile, 'link', quiet=TRUE)
gageIds <- paste0('g',formatC(rlLink,width=3,flag='0'))
newCopyId <- 'gagesAllLinks'
newFile <- AddRouteLinkGage(rlFile, gageIds, rlLink, new=newCopyId, overwrite=TRUE)
print(newFile)
## [1] "/Users/jamesmcc/WRF_Hydro/DOMAIN_library/Boulder_Creek_100m_1km_2sqkm_full_2015_09_03/Route_Link.gagesAllLinks.nc"
ncdump(newFile)
## File: /Users/jamesmcc/WRF_Hydro/DOMAIN_library/Boulder_Creek_100m_1km_2sqkm_full_2015_09_03/Route_Link.gagesAllLinks.nc
## ( NC_FORMAT_NETCDF4 ):
## dimensions (2):
## linkDim = 403 ;
## IDLength = 15 ;
## variables (26):
## integer IDLength(IDLength) ;
## IDLength:long_name = "IDLength"
## int link(linkDim) ;
## link:long_name = "Link ID" ;
## int from(linkDim) ;
## from:long_name = "From Link ID" ;
## from:coordinates = "time lat lon" ;
## int to(linkDim) ;
## to:long_name = "To Link ID" ;
## to:coordinates = "time lat lon" ;
## float lon(linkDim) ;
## lon:long_name = "longitude of the start node" ;
## lon:units = "degrees_east" ;
## lon:standard_name = "longitude" ;
## float lat(linkDim) ;
## lat:long_name = "latitude of the start node" ;
## lat:units = "degrees_north" ;
## lat:standard_name = "latitude" ;
## float alt(linkDim) ;
## alt:long_name = "Elevation in meters at start node" ;
## alt:standard_name = "height" ;
## alt:units = "m" ;
## alt:positive = "up" ;
## alt:axis = "Z" ;
## int type(linkDim) ;
## type:long_name = "Link type" ;
## type:coordinates = "time lat lon" ;
## int order(linkDim) ;
## order:long_name = "Stream order (Strahler)" ;
## order:coordinates = "time lat lon" ;
## float Qi(linkDim) ;
## Qi:long_name = "Initial flow in link (CMS)" ;
## Qi:coordinates = "time lat lon" ;
## float MusK(linkDim) ;
## MusK:long_name = "Muskingum routing time (s)" ;
## MusK:coordinates = "time lat lon" ;
## float MusX(linkDim) ;
## MusX:long_name = "Muskingum weighting coefficient" ;
## MusX:coordinates = "time lat lon" ;
## float Length(linkDim) ;
## Length:long_name = "Stream length (m)" ;
## Length:coordinates = "time lat lon" ;
## float n(linkDim) ;
## n:long_name = "Manning's roughness" ;
## n:coordinates = "time lat lon" ;
## float So(linkDim) ;
## So:long_name = "Slope (%; drop/length)" ;
## So:coordinates = "time lat lon" ;
## float ChSlp(linkDim) ;
## ChSlp:long_name = "Channel side slope (%; drop/length)" ;
## ChSlp:coordinates = "time lat lon" ;
## float BtmWdth(linkDim) ;
## BtmWdth:long_name = "Bottom width of channel" ;
## BtmWdth:coordinates = "time lat lon" ;
## float LkHZArea(linkDim) ;
## LkHZArea:long_name = "Maximum lake area (sq. m)" ;
## LkHZArea:coordinates = "time lat lon" ;
## float LkMxH(linkDim) ;
## LkMxH:long_name = "Maximum lake elevation (m)" ;
## LkMxH:coordinates = "time lat lon" ;
## float WeirC(linkDim) ;
## WeirC:long_name = "Weir coefficient" ;
## WeirC:coordinates = "time lat lon" ;
## float WeirL(linkDim) ;
## WeirL:long_name = "Weir length (m)" ;
## WeirL:coordinates = "time lat lon" ;
## float OrificeC(linkDim) ;
## OrificeC:long_name = "Orifice coefficient" ;
## OrificeC:coordinates = "time lat lon" ;
## float OrificeA(linkDim) ;
## OrificeA:long_name = "Orifice croess-sectional area (sq. m)" ;
## OrificeA:coordinates = "time lat lon" ;
## float OrificeE(linkDim) ;
## OrificeE:long_name = "Orifice elevation (m)" ;
## OrificeE:coordinates = "time lat lon" ;
## float time() ;
## time:standard_name = "time" ;
## time:long_name = "time of measurement" ;
## time:units = "days since 2000-01-01 00:00:00" ;
## char gages(linkDim,IDLength) ;
## gages:units = "usgsId" ;
## gages:_FillValue = " " ;
##
## // global attributes (2):
## :featureType = "point"
## :history = "Created Thu Sep 03 09:30:05 2015"
print(ncdump(newFile, 'gages', quiet=TRUE))
## [1] " g255" " g266" " g213"
## [4] " g253" " g277" " g208"
## [7] " g259" " g240" " g241"
## [10] " g257" " g217" " g216"
## [13] " g238" " g228" " g247"
## [16] " g229" " g231" " g258"
## [19] " g069" " g230" " g061"
## [22] " g244" " g063" " g222"
## [25] " g263" " g262" " g066"
## [28] " g256" " g252" " g115"
## [31] " g141" " g140" " g225"
## [34] " g103" " g067" " g065"
## [37] " g091" " g090" " g260"
## [40] " g187" " g227" " g146"
## [43] " g112" " g210" " g264"
## [46] " g133" " g132" " g163"
## [49] " g226" " g128" " g221"
## [52] " g092" " g215" " g110"
## [55] " g177" " g142" " g139"
## [58] " g234" " g199" " g164"
## [61] " g106" " g123" " g169"
## [64] " g118" " g209" " g170"
## [67] " g105" " g200" " g197"
## [70] " g126" " g223" " g186"
## [73] " g153" " g120" " g147"
## [76] " g175" " g107" " g179"
## [79] " g224" " g096" " g124"
## [82] " g154" " g089" " g151"
## [85] " g211" " g176" " g203"
## [88] " g138" " g167" " g099"
## [91] " g125" " g155" " g152"
## [94] " g150" " g181" " g207"
## [97] " g161" " g097" " g190"
## [100] " g058" " g280" " g278"
## [103] " g050" " g144" " g174"
## [106] " g098" " g188" " g059"
## [109] " g282" " g281" " g055"
## [112] " g113" " g204" " g104"
## [115] " g166" " g165" " g189"
## [118] " g121" " g056" " g052"
## [121] " g178" " g172" " g171"
## [124] " g300" " g295" " g101"
## [127] " g191" " g094" " g093"
## [130] " g316" " g122" " g057"
## [133] " g180" " g049" " g303"
## [136] " g302" " g325" " g095"
## [139] " g317" " g314" " g343"
## [142] " g053" " g212" " g040"
## [145] " g341" " g100" " g194"
## [148] " g087" " g082" " g332"
## [151] " g330" " g198" " g195"
## [154] " g319" " g086" " g083"
## [157] " g235" " g201" " g326"
## [160] " g318" " g085" " g205"
## [163] " g108" " g196" " g320"
## [166] " g088" " g310" " g116"
## [169] " g145" " g304" " g109"
## [172] " g074" " g232" " g068"
## [175] " g192" " g184" " g183"
## [178] " g084" " g305" " g296"
## [181] " g071" " g102" " g070"
## [184] " g292" " g288" " g315"
## [187] " g313" " g290" " g182"
## [190] " g080" " g078" " g173"
## [193] " g299" " g297" " g168"
## [196] " g162" " g289" " g157"
## [199] " g156" " g079" " g309"
## [202] " g308" " g143" " g075"
## [205] " g298" " g136" " g294"
## [208] " g285" " g248" " g119"
## [211] " g117" " g076" " g137"
## [214] " g261" " g291" " g322"
## [217] " g321" " g287" " g283"
## [220] " g249" " g329" " g242"
## [223] " g336" " g111" " g073"
## [226] " g038" " g324" " g323"
## [229] " g032" " g286" " g251"
## [232] " g345" " g254" " g340"
## [235] " g051" " g338" " g337"
## [238] " g062" " g339" " g328"
## [241] " g327" " g293" " g352"
## [244] " g348" " g347" " g346"
## [247] " g361" " g054" " g021"
## [250] " g243" " g018" " g047"
## [253] " g334" " g333" " g265"
## [256] " g335" " g359" " g006"
## [259] " g356" " g355" " g033"
## [262] " g127" " g378" " g376"
## [265] " g369" " g239" " g366"
## [268] " g365" " g362" " g360"
## [271] " g358" " g357" " g349"
## [274] " g214" " g077" " g012"
## [277] " g311" " g048" " g301"
## [280] " g039" " g034" " g030"
## [283] " g025" " g276" " g019"
## [286] " g273" " g268" " g267"
## [289] " g131" " g379" " g377"
## [292] " g375" " g374" " g245"
## [295] " g371" " g114" " g368"
## [298] " g367" " g237" " g364"
## [301] " g363" " g354" " g353"
## [304] " g350" " g284" " g014"
## [307] " g081" " g193" " g031"
## [310] " g312" " g046" " g041"
## [313] " g036" " g035" " g159"
## [316] " g028" " g027" " g023"
## [319] " g022" " g020" " g275"
## [322] " g015" " g270" " g269"
## [325] " g388" " g382" " g381"
## [328] " g380" " g373" " g372"
## [331] " g236" " g233" " g016"
## [334] " g351" " g220" " g218"
## [337] " g342" " g386" " g064"
## [340] " g060" " g307" " g306"
## [343] " g045" " g042" " g037"
## [346] " g160" " g389" " g029"
## [349] " g026" " g149" " g274"
## [352] " g017" " g272" " g013"
## [355] " g008" " g135" " g005"
## [358] " g387" " g130" " g385"
## [361] " g384" " g383" " g250"
## [364] " g246" " g370" " g391"
## [367] " g219" " g344" " g399"
## [370] " g206" " g331" " g202"
## [373] " g072" " g395" " g394"
## [376] " g185" " g393" " g044"
## [379] " g043" " g129" " g390"
## [382] " g158" " g024" " g279"
## [385] " g148" " g403" " g402"
## [388] " g401" " g400" " g271"
## [391] " g398" " g397" " g396"
## [394] " g011" " g010" " g009"
## [397] " g392" " g007" " g134"
## [400] " g004" " g003" " g002"
## [403] " g001"
Run the model with doubled precip
Set the symlink to the new Route_Link.nc file.
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek/DOMAIN> l
lrwxrwxrwx 2 jamesmcc rap 93 Sep 9 09:42 Route_Link.nc ->
../../DOMAIN_library/Boulder_Creek_100m_1km_2sqkm_full_2015_09_03/Route_Link.gagesAllLinks.nc
Set up several model parameters in the namelists.
namelist.hrldas:
START_YEAR = 2012
START_MONTH = 08
START_DAY = 01
START_HOUR = 00
START_MIN = 00
KHOUR = 48
FORCING_TIMESTEP = 3600
NOAH_TIMESTEP = 900
OUTPUT_TIMESTEP = 7200
FORC_TYP = 4
hydro.namelist:
out_dt = 15 ! minutes
DTRT_CH = 60
DTRT_TER = 6
SUBRTSWCRT = 1
OVRTSWCRT = 1
CHANRTSWCRT = 1
channel_option = 2
route_link_f = "DOMAIN/Route_Link.nc"
Run without nudging and using precip doubling. The wrf_hydro.noNudging_doublePrecip.exe binary was built with the following environment variables: PRECIP_DOUBLE=1 WRF_HYDRO_NUDGING=0 using the daBranch of wrf_hydro_model.
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> mpirun -np 1 ./wrf_hydro.noNudging_doublePrecip.exe
. . .
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> mkdir run.pmo
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> ./cleanup.sh run.pmo
Create timeslice files
chanObsFiles <- list.files('~/WRF_Hydro/Col_Bldr_Creek/run.pmo/', pattern = 'CHANOBS_DOMAIN', full.names=TRUE)
sliceResolutionMin <- 15
outputDir <- '~/WRF_Hydro/Col_Bldr_Creek/run.pmo/timeSlicePMO/'
dir.create(outputDir)
sliceFiles <- ChanObsToTimeSlice(chanObsFiles, sliceResolutionMin, outputDir)
Put these in place
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> rm nudgingTimeSliceObs
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> ln -sf run.pmo/timeSlicePMO/ nudgingTimeSliceObs
Create the nudging parameter file
min(ncdump(rlFile, 'Length', quiet=TRUE)) ## use R less than this /2.
## [1] 50
gageParams <- data.frame( gageId = paste0('g',formatC(1:403,width=3,flag='0')), stringsAsFactors=FALSE)
gageParams$R=24
gageParams$G=1
gageParams$tau=20
MkNudgingParams(gageId=gageParams$gageId, R=gageParams$R,
G=gageParams$G, tau=gageParams$tau,
outFile='~/WRF_Hydro/Col_Bldr_Creek/run.pmo/nudgingParams.PMOallGages.nc',
overwrite=TRUE)
## File: ~/WRF_Hydro/Col_Bldr_Creek/run.pmo/nudgingParams.PMOallGages.nc
## ( NC_FORMAT_NETCDF4 ):
## dimensions (2):
## stationIdStrLen = 15 ;
## stationIdInd = UNLIMITED ; // (403 currently)
## variables (4):
## char stationId(stationIdInd,stationIdStrLen) (Chunking: [15,1]) ;
## stationId:units = "-" ;
## stationId:long_name = "USGS station identifer" ;
## float R(stationIdInd) (Chunking: [1]) ;
## R:units = "meters" ;
## R:long_name = "Radius of influence in meters" ;
## float G(stationIdInd) (Chunking: [1]) ;
## G:units = "-" ;
## G:long_name = "Amplitude of nudging" ;
## float tau(stationIdInd) (Chunking: [1]) ;
## tau:units = "minutes" ;
## tau:long_name = "Time tapering parameter half window size in minutes" ;
Then put the param file in place.
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> ln -s run.pmo/nudgingParams.PMOallGages.nc nudgingParams.nc
Assimilate the observations using regular forcing
The wrf_hydro.nudging.exe binary was built with the following environment variables: PRECIP_DOUBLE=0 and WRF_HYDRO_NUDGING=1 using the daBranch of wrf_hydro_model.
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> mpirun -np 1 ./wrf_hydro.nudging.exe
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> mkdir run.pmoNudge
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> ./cleanup.sh run.pmoNudge
Run the model without nudging using regular forcingings.
The wrf_hydro.noNudging.exe binary was built with the following environment variables: PRECIP_DOUBLE=0 and WRF_HYDRO_NUDGING=0 using the daBranch of wrf_hydro_model.
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> mpirun -np 1 ./wrf_hydro.noNudging.exe
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> mkdir run.pmoNoNudge
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> ./cleanup.sh run.pmoNoNudge
Visualize
77779 frxst_pts_out.txt ## = 403stns*(48hours*4output/hr+1output@0)
frxstNoNudge <- ReadFrxstPts("~/WRF_Hydro/Col_Bldr_Creek/run.pmoNoNudge/frxst_pts_out.txt", stId='character')
frxstNudge <- ReadFrxstPts("~/WRF_Hydro/Col_Bldr_Creek/run.pmoNudge/frxst_pts_out.txt", stId='character')
## make sure the abscissae columns are equal before taking the differences
if(all(plyr::laply(NamedList(names(frxst)), function(colName) all(frxst[[colName]]==frxstNoNudge[[colName]]))[-6:-7]))
{
forcErr <- frxst
forcErr$`q err obs-model (cms)` <- frxst$q_cms - frxstNoNudge$q_cms
forcErr$`q err obs-model (cfs)` <- frxst$q_cfs - frxstNoNudge$q_cfs
forcErr$q_cms <- forcErr$q_cfs <- NULL
forcErr$error <- 'forcing'
} else warning('Abscissae do not line up with frxstNoNudge, please investigate', immediate.=TRUE)
if(all(plyr::laply(NamedList(names(frxst)), function(colName) all(frxst[[colName]]==frxstNudge[[colName]]))[-6:-7]))
{
nudgeErr <- frxst
nudgeErr$`q err obs-model (cms)` <- frxst$q_cms - frxstNudge$q_cms
nudgeErr$`q err obs-model (cfs)` <- frxst$q_cfs - frxstNudge$q_cfs
nudgeErr$q_cms <- nudgeErr$q_cfs <- NULL
nudgeErr$error <- 'nudging'
} else warning('Abscissae do not line up with frxstNudge, please investigate', immediate.=TRUE)
err <- rbind(forcErr, nudgeErr)
ggplot2::ggplot(subset(err, POSIXct < as.POSIXct('2012-08-01 06:00:00', tz='UTC')), ## all the action is early on
ggplot2::aes(x=POSIXct, y=`q err obs-model (cms)`, color=st_id)) +
ggplot2::geom_line() +
ggplot2::scale_color_discrete(guide='none') +
ggplot2::facet_wrap(~error, ncol=1) +
ggplot2::theme_bw()
